Setup

Set up environment and load in data

Adding in parameters from the two systems model for all subjects (i.e. not choosing the best fitting one per subject)

source(paste0(helpers_path,'ddModels/sim_task.R')) # do this first not to mess with cluster setup?
source(paste0(helpers_path,'ddModels/sim_sanity_checks.R'))
source(paste0(helpers_path,'rlModels/fit_rl_hierarchical.R'))
## [1] "Beginning reading data in..."
## [1] "Done!"
source(paste0(helpers_path,'add_inferred_pars.R'))
clean_beh_data = add_inferred_pars(clean_beh_data, par_ests, model_name="original")

Create empty list that will store the trial simulators for the forthcoming models.

sim_trial_list = list()

True data

Filter data for a few subjects with a range of learning rates (i.e. variance in QValue differences) to simulate RTs using DDM. Rename columns to work with task simulation function.

sub_data = clean_beh_data %>%
  filter(subnum  %in% c("01", "03", "05","07", "09", "11", "13", "15", "17", "19")) %>%
  select(leftQValue, rightQValue, leftLotteryEV, rightLotteryEV, probFractalDraw, reactionTime, choiceLeft, subnum) %>%
  rename(EVLeft = leftLotteryEV, EVRight = rightLotteryEV, QVLeft = leftQValue, QVRight = rightQValue)

Sanity checks in the sampled data

sim_sanity_checks(sub_data %>%
                    select(-subnum) %>%
                    mutate(choice = ifelse(choiceLeft == 1, "left", "right")), 
                  compare_rts = FALSE)
## [1] "Proportion of time out trials if no decision made: 0"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Individual subject plots

For the above stylized RT and choice plots

sub_data %>%
  # select(-subnum) %>%
  mutate(choice = ifelse(choiceLeft == 1, "left", "right"))%>%
  drop_na()%>%
  mutate(probFractalDraw = as.factor(probFractalDraw), 
         log_rt = log(reactionTime)) %>%
  group_by(probFractalDraw, subnum) %>%
  summarise(mean_log_rt = mean(log_rt),
            sem_log_rt = sem(log_rt), .groups="keep") %>%
  ggplot(aes(probFractalDraw, mean_log_rt))+
  geom_point()+
  geom_errorbar(aes(ymin = mean_log_rt - sem_log_rt, ymax = mean_log_rt + sem_log_rt), width=.2)+
  facet_wrap(~subnum, scales = 'free')

sub_data %>%
  mutate(choice = ifelse(choiceLeft == 1, "left", "right"))%>%
  select(EVLeft, EVRight, QVLeft, QVRight, probFractalDraw, choice, reactionTime, subnum) %>%
  mutate(probFractalDraw = as.factor(probFractalDraw),
         choiceLeft = ifelse(choice == "left", 1, ifelse(choice=="right", 0, NA)),
         EVDiff = EVLeft - EVRight, 
         QVDiff = QVLeft - QVRight) %>%
  nest(data = -probFractalDraw, -subnum) %>% 
  mutate(
    fit = map(data, ~ glm(choiceLeft ~ scale(EVDiff) + scale(QVDiff), data = .x, family=binomial(link="logit"))),
    tidied = map(fit, tidy)
  ) %>% 
  unnest(tidied) %>%
  filter(term != "(Intercept)") %>%
  select(subnum, probFractalDraw, term, estimate, std.error) %>%
  ggplot(aes(probFractalDraw, estimate, col=term, group=term))+
  geom_point()+
  geom_line()+
  geom_errorbar(aes(ymin = estimate - std.error, ymax = estimate +std.error), width=0.2)+
  geom_hline(aes(yintercept=0), linetype="dashed")+
  scale_color_manual(values = cbbPalette[2:1])+
  theme(legend.position = "bottom")+
  labs(color="", title="Relevant attribute effect on choice")+
  facet_wrap(~subnum, scales="free")+
  ylim(-2.5, 10)

Model 1: Simplest

  • Integration begins at stim presentation for all conditions
  • Drift rate is proportional to the bundle value difference
  • Bundle values are computed as sums of QV and EV weighted by their relevance (no distortion of probability)
source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model1.R'))
sim_trial_list[['model1']] = sim_trial

Check sim task

m1 = sim_task(sub_data, model_name="model1", d=0.04, sigma = 0.02)
sim_sanity_checks(m1)
## [1] "Proportion of time out trials if no decision made: 0.215"
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: Removed 639 rows containing non-finite values (stat_bin).

## Warning: glm.fit: algorithm did not converge

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning: glm.fit: algorithm did not converge
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

## Warning: Removed 18 rows containing missing values (geom_point).
## Warning: Removed 18 row(s) containing missing values (geom_path).

m1_2 = sim_task(sub_data, model_name="model1", d=0.04, sigma = 0.02, epsilon = 0.03)
sim_sanity_checks(m1_2)

Spread vs mean

Model 2: Early integration

  • Early integrator for pFractalDraw == 1 only.
source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model2.R'))
sim_trial_list[['model2']] = sim_trial

Check sim task

m2 = sim_task(sub_data, model_name = "model2", d=0.05, sigma = 0.02)
sim_sanity_checks(m2)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.054"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.165"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

## Warning: Removed 2 rows containing missing values (geom_point).

When does integration reach a bound before the stim presentation?

with(m2, table(decPreStim, probFractalDraw))
##           probFractalDraw
## decPreStim   0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1
##          0 244 249 247 294 297 297 297 299 248 250  92
##          1   2   0   1   0   0   0   2   0   1   0 154

RT distributions depending on whether a decision was reached before stimulus presentation.

m2 %>%
  select(EVLeft, EVRight, QVLeft, QVRight, probFractalDraw, choice, reactionTime, decPreStim) %>%
  mutate(data_type = "sim") %>%
  rbind(sub_data %>%
          mutate(choice = ifelse(choiceLeft == 1, "left", "right"),
                 data_type = "true",
                 decPreStim = "trueData") %>%
          select(-subnum, -choiceLeft)) %>%
  mutate(probFractalDraw = as.factor(probFractalDraw)) %>%
  ggplot(aes(reactionTime, fill=decPreStim)) +
  geom_histogram(position="identity", bins=30, alpha=.5) +
  facet_wrap(~probFractalDraw)+
  theme(legend.position = "bottom")

Are the number of timeout trials similar across conditions? Yes.

with(m2, table(timeOut, probFractalDraw))
##        probFractalDraw
## timeOut   0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9   1
##       0 192 210 208 256 253 252 241 254 200 200 218
##       1  54  39  40  38  44  45  58  45  49  50  28

Model 2a: asymmetric weighting

  • Single parameter theta controls how much lotteries should be overweighted and how much fractals should be underweighted. It is an initial effort to induce different drift rates for the attributes.
source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model2a.R'))
sim_trial_list[['model2a']] = sim_trial

Check sim task

m2a = sim_task(sub_data, model_name = "model2a", d=0.05, sigma = 0.02, theta = 0.025)
sim_sanity_checks(m2a, yrange_lim = 30)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.053"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.165"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

## Warning: Removed 2 rows containing missing values (geom_point).

Model 2b: Asymmetric prob distortion

  • Distort probFractalDraw when integrating info about fractals but no distortion for lotteries
  • Intended to capture the stepwise nature of the logit slopes for the QV difference but the linear nature for the EV difference
data.frame(pFrac = seq(0, 1, .1),
           delta = 3, gamma = 3) %>%
  mutate(distortedPFrac = exp( (-1) * delta * ((-1)*log(pFrac))^gamma) ) %>%
  ggplot(aes(pFrac, distortedPFrac))+
  geom_point()+
  geom_line()+
  geom_abline(aes(slope=1, intercept=0), linetype="dashed")+
  scale_y_continuous(breaks = seq(0, 1, .1))+
  scale_x_continuous(breaks = seq(0, 1, .1))

Overweight for larger pfrac

source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model2b.R'))
sim_trial_list[['model2b']] = sim_trial

Check sim task

Can we get 0 logit slopes for probFractalDraw < .5 trials for the QV difference with an assymmetric prob weighting curve as above? Yes.

m2b_1 = sim_task(sub_data, model_name = "model2b", d=0.04, sigma = 0.01, delta = 3, gamma = 3)
sim_sanity_checks(m2b_1, yrange_lim = 50)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.048"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.213"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

## Warning: Removed 1 rows containing missing values (geom_point).

Can we reduce timeouts with a barrier decay? Yes but not entirely.

m2b_2 = sim_task(sub_data, model_name = "model2b", d=0.03, sigma = 0.01, delta = 2, gamma = 2, barrierDecay = .007)
sim_sanity_checks(m2b_2, yrange_lim = 50)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.065"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.055"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Why are the logit values so high?

In the simulations choice depends only on these value differences. The model fits too well, there is very low residual deviance. Here’s the true data logit:

tmp = sub_data %>%
  filter(probFractalDraw == 1) %>%
  mutate(QVDiff = scale(QVLeft - QVRight),
         EVDiff = scale(EVLeft - EVRight))

summary(glm(choiceLeft ~ EVDiff + QVDiff, data = tmp,family=binomial(link="logit")))
## 
## Call:
## glm(formula = choiceLeft ~ EVDiff + QVDiff, family = binomial(link = "logit"), 
##     data = tmp)
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -2.3305  -0.8748  -0.2015   0.9356   2.3326  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -0.1419     0.1503  -0.944    0.345    
## EVDiff        0.1366     0.1489   0.917    0.359    
## QVDiff        1.5812     0.2354   6.718 1.84e-11 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 340.23  on 245  degrees of freedom
## Residual deviance: 262.26  on 243  degrees of freedom
## AIC: 268.26
## 
## Number of Fisher Scoring iterations: 5

And the logit for the simulation

tmp = m2b_2 %>%
  filter(probFractalDraw == 1) %>%
  mutate(QVDiff = scale(QVLeft - QVRight),
         EVDiff = scale(EVLeft - EVRight),
         choiceLeft = ifelse(choice == "left", 1, 0))

summary(glm(choiceLeft ~ EVDiff + QVDiff, data = tmp,family=binomial(link="logit")))
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## 
## Call:
## glm(formula = choiceLeft ~ EVDiff + QVDiff, family = binomial(link = "logit"), 
##     data = tmp)
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -1.70314  -0.00032   0.00000   0.00000   2.52333  
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)   
## (Intercept)  -5.1736     1.7955  -2.881  0.00396 **
## EVDiff       -0.3134     0.6089  -0.515  0.60678   
## QVDiff       35.4864    12.0201   2.952  0.00315 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 331.602  on 245  degrees of freedom
## Residual deviance:  24.371  on 243  degrees of freedom
## AIC: 30.371
## 
## Number of Fisher Scoring iterations: 12
  • Can we reduce the logit slopes by increasing epsilon? (SD of distribution the mean of drift is sampled from) Yes.

But does it also destroy the RT inverse U? Do the average RTs become more similar to each other for all intermediate probFractalDraw levels?

m2b_3 = sim_task(sub_data, model_name = "model2b", d=0.02, sigma = 0.007, delta = 3, gamma = 3, barrierDecay = 0.004, epsilon = 0.03, nonDecisionTime = 0)
sim_sanity_checks(m2b_3)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.345"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.008"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

  • What do the trials on the slow end of the bimodal distribution look like? Are they for particularly hard trials? Yes, for slow trials the difference between the options is always lower.

Should we add something to an effect of “if options are too close to each other choose…” randomly? the first one?

m2b_3 %>%
  mutate(absQVDiff= abs(QVLeft - QVRight),
         absEVDiff = abs(EVLeft - EVRight),
         slowTrial = factor(ifelse(reactionTime > 3, 1, 0))) %>%
  select(absQVDiff, absEVDiff, slowTrial) %>%
  gather(key, value, -slowTrial)%>%
  group_by(key, slowTrial) %>%
  summarise(.groups = 'keep',
            mean_diff = mean(value),
            sem_diff = sem(value)) %>%
  ggplot(aes(key, mean_diff, color=slowTrial))+
  geom_point(position=position_dodge(width = .5))+
  geom_errorbar(aes(ymin = mean_diff-sem_diff, ymax= mean_diff + sem_diff), width = .2, position=position_dodge(width = .5))+
  labs(x="")+
  theme(legend.position = "bottom")

Model 4a: 3 integrators

  • 3 integrators all with their own drift and noise rates
  • Arbitrator integrator starts biased towards lotteries
  • If pFractalDraw == 1 fractal integrator (not the arbitrator) starts biased towards the better fractal
  • Attribute relevance (pFractalDraw) affects attribute integrators directly, not the arbitrator
  • Arbitrator depends on the difference in absolute attribute integrator RDVs
source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model4a.R'))
sim_trial_list[['model4a']] = sim_trial

Check sim task

m4a_1 = sim_task(sub_data, model_name = "model4a", dLott=0.03, dFrac=0.04, dArb=0.04, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01)
sim_sanity_checks(m4a_1)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.029"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Reduce logit slopes and try to increase lottery bias to see if you would see a change in logit slopes at pFrac == 0.5

m4a_2 = sim_task(sub_data, model_name = "model4a", dLott=0.025, dFrac=0.04, dArb=0.045, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, epsilon = 0.003, lotteryBias = 0.35)
sim_sanity_checks(m4a_2)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.026"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Model 4b: 4a with asymmetric prob distortion

source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model4b.R'))
sim_trial_list[['model4b']] = sim_trial

Check sim task

m4b_1 = sim_task(sub_data, model_name = "model4b", dLott=0.03, dFrac=0.04, dArb=0.04, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, delta = 3, gamma = 3)
sim_sanity_checks(m4b_1)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.03"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Add some noise to decrease logit slopes

m4b_2 = sim_task(sub_data, model_name = "model4b", dLott=0.03, dFrac=0.04, dArb=0.04, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, delta = 3, gamma = 3, epsilon = 0.015)
sim_sanity_checks(m4b_2)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.034"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Add barrier decay to reduce time outs

m4b_3 = sim_task(sub_data, model_name = "model4b", dLott=0.025, dFrac=0.035, dArb=0.035, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, delta = 3, gamma = 1.5, epsilon = 0.015, barrierDecay = 0.003)
sim_sanity_checks(m4b_3)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.01"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Model 5a: 3 integrators

  • 3 integrators all with their own drift and noise rates
  • Arbitrator integrator starts biased towards lotteries
  • If pFractalDraw == 1 fractal integrator (not arbitrator) starts biased towards the better fractal
  • Attribute relevant (pFractalDraw) does not affects attribute integrators directly
  • Arbitrator integrator depends on the difference in absolute attribute integrator RDVs weighted by the attribute relevance
source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model5a.R'))
sim_trial_list[['model5a']] = sim_trial

Example visualization of integration

Check sim task

m5a_1 = sim_task(sub_data, model_name = "model5a", dLott=0.03, dFrac=0.04, dArb=0.04, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01)
sim_sanity_checks(m5a_1)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.037"
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

## Warning: Removed 2 rows containing missing values (geom_point).
## Warning: Removed 2 row(s) containing missing values (geom_path).

Slow down lottery integrator and lower the logit slopes

m5a_2 = sim_task(sub_data, model_name = "model5a", dLott=0.02, dFrac=0.04, dArb=0.04, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, epsilon = 0.025)
sim_sanity_checks(m5a_2)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.033"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

RT distributions based on arbitrator used for the decision

m5a_2 %>%
  select(EVLeft, EVRight, QVLeft, QVRight, probFractalDraw, choice, reactionTime, arbitrator) %>%
  mutate(data_type = "sim") %>%
  rbind(sub_data %>%
          mutate(choice = ifelse(choiceLeft == 1, "left", "right"),
                 data_type = "true",
                 arbitrator= "true") %>%
          select(-subnum, -choiceLeft)) %>%
  mutate(probFractalDraw = as.factor(probFractalDraw)) %>%
  ggplot(aes(reactionTime, fill=as.factor(arbitrator))) +
  geom_histogram(position="identity", bins=30, alpha=.5) +
  labs(title="RT long tail?", fill="")+
  facet_wrap(~probFractalDraw)+
  theme(legend.position = "bottom")

Model 5b: 5a with asymmetric prob distortion

source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model5b.R'))
sim_trial_list[['model5b']] = sim_trial

Check sim task

m5b_1 = sim_task(sub_data, model_name = "model5b", dLott=0.02, dFrac=0.04, dArb=0.035, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, delta = 3, gamma = 3, epsilon = .02)
sim_sanity_checks(m5b_1)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.031"

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Model 6a: 5a + early integration

  • 3 integrators all with their own drift and noise rates
  • Arbitrator integrator starts from lottery bias
  • If pFractalDraw == 1 fractal and arbitrator integrator starts early. Arbitrator moves only towards fractal.
  • Attribute relevance (pFractalDraw) does not affect attribute integrators directly
  • Arbitrator integrator depends on the difference in absolute attribute integrator RDVs weighted by attribute relevance
source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model6a.R'))
sim_trial_list[['model6a']] = sim_trial

Check sim task

m6a = sim_task(sub_data, model_name = "model6a", dLott=0.015, dFrac=0.02, dArb=0.05, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, lotteryBias = .25, epsilon = 0.003)
sim_sanity_checks(m6a)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.079"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.049"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Model 6b: 6a with asymmetric prob distortion

source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model6b.R'))
sim_trial_list[['model6b']] = sim_trial

Check sim task

m6b_1 = sim_task(sub_data, model_name = "model6b", dLott=0.02, dFrac=0.04, dArb=0.04, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, delta = 3, gamma = 3)
sim_sanity_checks(m6b_1)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.078"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.034"

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred

## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Model 6c: 4b with early integration

source(paste0(helpers_path, 'ddModels/r_ddm_models/ddm_model6c.R'))
sim_trial_list[['model6c']] = sim_trial

Check sim task

m6c_1 = sim_task(sub_data, model_name = "model6c", dLott=0.03, dFrac=0.04, dArb=0.04, sigmaLott = 0.03, sigmaFrac = 0.03, sigmaArb = 0.01, delta = 3, gamma = 3)
sim_sanity_checks(m6c_1)
## [1] "Proportion of time out trials if no decision made: 0"
## [1] "Proportion of trials with decision before stim: 0.075"
## [1] "Proportion of time out trials with decision of closest boundary sampled RTs: 0.029"

## Warning in if (large_range) {: the condition has length > 1 and only the first
## element will be used

Stop cluster

parallel::stopCluster(cl = my.sim.cluster)
# rm(my.cluster)

Trade-offs across models

  • Is the relationship between pars that optim rt vs choice the same across models?

Older models

Model 3: Model 2 + relative preference weighting

  • Early integration for pFractalDraw == 1 only
  • Relative preference weighted value difference for all other conditions

Model 4, 5, 6, 7

  • Same as the versions with the ‘a’ suffix, except lotteries and fractals have the same drift rate and noise level

Model 7a: Switching between integrators

  • 3 integrators all with their own drift and noise rates
  • Arbitrator integrator starts from lottery bias unless pFractalDraw == 1. Then it starts from 0.
  • If pFractalDraw == 1 fractal and arbitrator integrator starts early. Arbitrator moves only towards fractal.
  • Attribute relevance (pFractalDraw) does not affect attribute integrators directly
  • Arbitrator integrator depends on the difference in absolute attribute integrator RDVs weighted by attribute relevance
  • When stimulus is presented arbitrator moves only based on lotteryRDV for a short period

Questions

  • True data logit slopes for QV difference never exceeds the slopes of lotteries. Why do I get this pattern in simulations?
  • Why is the slowest condition pFractalDraw == .6 for simulations but pFractalDraw == .5 in true data?
  • Where should NDT come into play for models with early integration?
  • Change grid search optimization for RT to KL divergence?
  • Incorporating EV and QV computation instead of using them as input. If these models are capturing what is happening after the prob fractal draw screen until a decision is made then QV can remain as an input (maybe with an NDT for recall) but EV still can only be computed during the stimulus screen.